BB 16S MB Final Analyses

Load required packages.

require(phyloseq)

Loading required package: phyloseq

 require(parallelDist)

Loading required package: parallelDist
Warning: package ‘parallelDist’ was built under R version 4.1.2

require(ggplot2)

Loading required package: ggplot2
Warning: package ‘ggplot2’ was built under R version 4.1.2

  require(randomForest)

Loading required package: randomForest
Warning: package ‘randomForest’ was built under R version 4.1.2
randomForest 4.7-1.1
Type rfNews() to see new features/changes/bug fixes.

Attaching package: ‘randomForest’

The following object is masked from ‘package:ggplot2’:

    margin

require(Boruta)

Loading required package: Boruta
Warning: package ‘Boruta’ was built under R version 4.1.2

  require(phyloseq)
require(stats)
  require(dplyr)

Loading required package: dplyr
Warning: package ‘dplyr’ was built under R version 4.1.2

Attaching package: ‘dplyr’

The following object is masked from ‘package:randomForest’:

    combine

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

require(readr)

Loading required package: readr
Warning: package ‘readr’ was built under R version 4.1.2

  require(pheatmap)

Loading required package: pheatmap

require(tidyverse)

Loading required package: tidyverse
Warning: package ‘tidyverse’ was built under R version 4.1.2
Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      

── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.2 ──
✔ tibble  3.1.8     ✔ stringr 1.5.0
✔ tidyr   1.2.1     ✔ forcats 0.5.2
✔ purrr   0.3.5     

Warning: package ‘tibble’ was built under R version 4.1.2
Warning: package ‘tidyr’ was built under R version 4.1.2
Warning: package ‘purrr’ was built under R version 4.1.2
Warning: package ‘stringr’ was built under R version 4.1.2
Warning: package ‘forcats’ was built under R version 4.1.2

── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::combine()       masks randomForest::combine()
✖ dplyr::filter()        masks stats::filter()
✖ dplyr::lag()           masks stats::lag()
✖ randomForest::margin() masks ggplot2::margin()

  require(micrUBIfuns)

Loading required package: micrUBIfuns
Loading required package: data.table
Warning: package ‘data.table’ was built under R version 4.1.2
data.table 1.14.6 using 1 threads (see ?getDTthreads).  Latest news: r-datatable.com
**********
This installation of data.table has not detected OpenMP support. It should still work but in single-threaded mode.
This is a Mac. Please read https://mac.r-project.org/openmp/. Please engage with Apple and ask them for support. Check r-datatable.com for updates, and our Mac instructions here: https://github.com/Rdatatable/data.table/wiki/Installation. After several years of many reports of installation problems on Mac, it's time to gingerly point out that there have been no similar problems on Windows or Linux.
**********

Attaching package: ‘data.table’

The following object is masked from ‘package:purrr’:

    transpose

The following objects are masked from ‘package:dplyr’:

    between, first, last

Loading required package: ggsci
Loading required package: ComplexHeatmap
Loading required package: grid
========================================
ComplexHeatmap version 2.8.0
Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
Github page: https://github.com/jokergoo/ComplexHeatmap
Documentation: http://jokergoo.github.io/ComplexHeatmap-reference

If you use it in published research, please cite:
Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
  genomic data. Bioinformatics 2016.

The new InteractiveComplexHeatmap package can directly export static 
complex heatmaps into an interactive Shiny app with zero effort. Have a try!

This message can be suppressed by:
  suppressPackageStartupMessages(library(ComplexHeatmap))
========================================
! pheatmap() has been masked by ComplexHeatmap::pheatmap(). Most of the arguments
   in the original pheatmap() are identically supported in the new function. You 
   can still use the original function by explicitly calling pheatmap::pheatmap().


Attaching package: ‘ComplexHeatmap’

The following object is masked from ‘package:pheatmap’:

    pheatmap

Loading required package: vegan
Warning: package ‘vegan’ was built under R version 4.1.2
Loading required package: permute
Warning: package ‘permute’ was built under R version 4.1.2
Loading required package: lattice
This is vegan 2.6-4

require(circlize)

Loading required package: circlize
Warning: package ‘circlize’ was built under R version 4.1.2
========================================
circlize version 0.4.15
CRAN page: https://cran.r-project.org/package=circlize
Github page: https://github.com/jokergoo/circlize
Documentation: https://jokergoo.github.io/circlize_book/book/

If you use it in published research, please cite:
Gu, Z. circlize implements and enhances circular visualization
  in R. Bioinformatics 2014.

This message can be suppressed by:
  suppressPackageStartupMessages(library(circlize))
========================================

  require(phylosmith)

Loading required package: phylosmith

require(phangorn)

Loading required package: phangorn
Warning: package ‘phangorn’ was built under R version 4.1.2
Loading required package: ape
Warning: package ‘ape’ was built under R version 4.1.2

Attaching package: ‘ape’

The following object is masked from ‘package:circlize’:

    degree


Attaching package: ‘phangorn’

The following objects are masked from ‘package:vegan’:

    diversity, treedist

  require(gridExtra)

Loading required package: gridExtra

Attaching package: ‘gridExtra’

The following object is masked from ‘package:dplyr’:

    combine

The following object is masked from ‘package:randomForest’:

    combine

require(pairwiseAdonis)

Loading required package: pairwiseAdonis
Loading required package: cluster
Warning: package ‘cluster’ was built under R version 4.1.2

  require(car)

Loading required package: car
Warning: package ‘car’ was built under R version 4.1.2
Loading required package: carData
Warning: package ‘carData’ was built under R version 4.1.2

Attaching package: ‘car’

The following object is masked from ‘package:purrr’:

    some

The following object is masked from ‘package:dplyr’:

    recode

Load necessary files from Dada2 processing script.

load("/Users/gordoncuster/Desktop/Git_Projects/BeeBread/Data/16S/DADA2outputs/dada2_outputs_16S_BeeBread_wtree.RData")

Read in metadata

metadata<-read.csv("/Users/gordoncuster/Desktop/Git_Projects/BeeBread/Data/BB_sample-metadata.txt", sep = "\t")
rownames(metadata)<-metadata$X
metadata$Region_Site<-paste( metadata$Region, metadata$Site, sep = "")
metadata$Treatment_Site<-paste( metadata$Treatment, metadata$Site, sep = "")

Convert files to correct format for phyloseq object and then rename samples in OTU table to match the format of the metadata file. Create phyloseq object and root phylogenetic tree for usage with distance metrics (e.g., Unifrac).

#assign object type to merge into phyloseq object
md<-sample_data(metadata)
otu<-otu_table(seqtab.nochim, taxa_are_rows = F)
tax_tab<-tax_table(taxa)
#check sample names of both otu table and metadata
#remove .fastq.gz from otu table names
sample_names(otu)<-str_split(sample_names(otu), pattern = ".fastq.gz", simplify = T) [,1]
#since they match, you can create the phyloseq object
bb16S_orig<-phyloseq(md, otu, tax_tab, fitGTR$tree)

#root tree for distance metrics like unifrac
set.seed(11)
phy_tree(bb16S_orig)<-root(phy_tree(bb16S_orig), sample(taxa_names(bb16S_orig), 1), resolve.root = TRUE)
is.rooted(phy_tree(bb16S_orig))

[1] TRUE

#Pre-processing. Rename OTU IDS to something more manageable. We save the ASCV sequences in case we need them later (e.g., BLASTN searches). Create Bacterial and Archaeal phyloseq objects, but first we remove Chloroplasts.

#Extract original IDs for future reference. 
seq_df_w_OTUID<-data.frame(OTUID = 1:ntaxa(bb16S_orig), Sequence = taxa_names(bb16S_orig))
taxa_names(bb16S_orig)<-paste("OTU_" , 1:ntaxa(bb16S_orig), sep = "")


#remove chloroplast
#978 chloroplast or mitochondrial taxa
bb16S_orig_nc<-subset_taxa(bb16S_orig, Family!= "Mitochondria" | is.na(Family) & Order!="Chloroplast" | is.na(Order))


#1337 bacterial taxa - 9 non-bacteria removed. 
bb16S_bac<-subset_taxa(bb16S_orig_nc, Kingdom == "Bacteria")
#9 archaea
bb16S_arch<-subset_taxa(bb16S_orig_nc, Kingdom == "Archaea")

Remove control samples as there is only as single replicate of each.

ps_wo_control<-subset_samples(bb16S_bac, Treatment != "none")

Examine rarefaction curves. Plateau in most samples under the 1.5k read mark.

rarecurve(otu_table(ps_wo_control), step=50, cex=0.5, main = "Rarefaction Curve", ylab="Richness", xlab="Read Depth")

Error in as(x, "matrix")[i, j, drop = FALSE] : 
  (subscript) logical subscript too long

set.seed(11)
sort(sample_sums(ps_wo_control))

  Dal_IPM1054_A   Dal_CON1058_A   Dal_CON1057_A   Spr_CON3017_A   Dal_IPM1055_A    Dal_CF1052_A    Spr_CF3016_A   Dal_IPM1055_B   Dal_CON1058_B   Sum_IPM4002_B 
            409             496             894            1188            1378            1469            1718            1962            2036            2037 
  Spr_CON3017_B    Sla_CF2066_B   Sla_IPM2071_B    Dal_CF1051_A    Dal_CF1051_B   Dal_IPM1054_B   Dal_CON1057_B    Sla_CF2065_B   Spr_CON3019_A   Sla_IPM2070_A 
           2254            2287            2288            2318            2359            2370            2396            2744            2835            2864 
Spr_CON3019_Be2    Spr_CF3013_A   Spr_IPM3023_B   Sla_CON2061_A   Spr_IPM3023_A   Sum_CON4011_B   Bel_IPM3042_A   Sla_IPM2070_B    Dal_CF1052_B   Sum_IPM4001_A 
           2886            3503            3513            3871            3898            4039            4196            4553            4710            4721 
   Sla_CF2066_A    Bel_CF3047_A   Sum_CON4011_A   Sum_IPM4001_B    Spr_CF3016_B    Sum_CF4006_B    Sum_CF4005_A    Bel_CF3045_B    Sla_CF2065_A   Sum_IPM4002_A 
           4749            4845            5502            5625            5693            5819            5935            5938            6648            6815 
  Spr_IPM3024_A    Ber_CF2006_B   Bel_CON3039_A   Sum_CON4010_B   Sla_CON2063_A   Ber_IPM2002_B   Spr_IPM3024_B    Spr_CF3013_B   Bel_CON3037_A    Sum_CF4006_A 
           6930            7047            7119            7147            7253            7462            7669            7888            8129            8276 
  Bel_IPM3042_B Ber_CON2011_Be2    Sum_CF4005_B    Bel_CF3045_A   Sum_CON4010_A   Ber_CON2010_B   Ber_CON2011_A   Bel_IPM3044_B   Sla_IPM2071_A   Bel_CON3039_B 
           9187            9315            9820           10444           10531           11063           11158           11247           12197           12295 
  Sla_CON2061_B    Bel_CF3047_B    Ber_CF2007_A   Bel_IPM3044_A    Ber_CF2007_B   Ber_IPM2002_A   Ber_CON2010_A   Bel_CON3037_B   Ber_IPM2001_A   Ber_IPM2001_B 
          13463           14099           15170           15601           17086           18639           19194           20875           21680           22834 
   Ber_CF2006_A   Sla_CON2063_B 
          33732           47127 

summary(sample_sums(ps_wo_control))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    409    2812    5877    8020   10466   47127 

sd(sample_sums(ps_wo_control))

[1] 7801.64

bb16S_bac_rarefy<-rarefy_even_depth(ps_wo_control, sample.size = 1100, rngseed = 14, trimOTUs = T)

`set.seed(14)` was used to initialize repeatable random subsampling.
Please record this for your records so others can reproduce.
Try `set.seed(14); .Random.seed` for the full vector
...
3 samples removedbecause they contained fewer reads than `sample.size`.
Up to first five removed samples are: 

Dal_CON1057_ADal_CON1058_ADal_IPM1054_A
...
486OTUs were removed because they are no longer 
present in any sample after random subsampling

...

#517 otus were removed during rarefaction. 

#To move forward, we can hellinger transform and go from there.
bb16S_bac_hellinger<-transform_sample_counts(ps_wo_control, function(x) sqrt(x / sum(x)))

#maximum liklihood point estimates
bb16S_bac_ML<- transform_sample_counts(ps_wo_control, function(x) x / sum(x))

#Alpha Diversity

sf2a <- plot_richness(bb16S_bac_rarefy, x = "Site", measures = c("Shannon", "Observed", "Chao1"), color = "Site") + geom_boxplot() + ggtitle("Sampling Location") + theme_classic()  + xlab("Sampling Location")  + theme(axis.text=element_text(size=12, angle = 45, hjust = 0.75),
        axis.title=element_text(size=14))

sf2b <- plot_richness(bb16S_bac_rarefy, x = "Site", measures = c("Shannon", "Observed", "Chao1"), color = "Treatment") + geom_boxplot() + ggtitle("Sampling Location by Treatment") + theme_classic() + ylab("") + xlab("Sampling Location") +theme(axis.text=element_text(size=12, angle = 45, hjust = 0.75),
        axis.title=element_text(size=14))

sf2c <-plot_richness(bb16S_bac_rarefy, x = "Site", measures = c("Shannon", "Observed", "Chao1"), color = "Timepoint") + geom_boxplot() + ggtitle("Sampling Location by Sammpling Time") + theme_classic() + xlab("Sampling Location")  +theme(axis.text=element_text(size=12, angle = 45, hjust = 0.75),
        axis.title=element_text(size=14))

sf2d <-plot_richness(bb16S_bac_rarefy, x = "Timepoint", measures = c("Shannon", "Observed", "Chao1"), color = "Treatment") + geom_boxplot() + ggtitle("Treatment by Sampling Time") + theme_classic() + ylab("") + xlab("Sampling Timepoint")  +theme(axis.text=element_text(size=12, angle = 45, hjust =0.75),
        axis.title=element_text(size=14))

sf2 <- grid.arrange(sf2a, sf2b, sf2c, sf2d)

plot(sf2)

Calculate alpha diversity metrics

#create alpha diversity table and prep data to include grouping columns
richness_BB16S_bac<-estimate_richness(bb16S_bac_rarefy, measures = c("Shannon", "Observed", "Chao1"))
#add metadata
richness_BB16S_bac$Site<-sample_data(bb16S_bac_rarefy)$Site
richness_BB16S_bac$Timepoint<-sample_data(bb16S_bac_rarefy)$Timepoint
richness_BB16S_bac$Treatment<-sample_data(bb16S_bac_rarefy)$Treatment
richness_BB16S_bac$Treatment_Timepoint<-sample_data(bb16S_bac_rarefy)$Treatment_Timepoint

Tests of Shannon diversity

#ANOVA and Tukey's PW comparisons
#Run a full 3 way ANOVA with interactions
mod<-aov(Shannon ~ Site * Treatment * Timepoint, data = richness_BB16S_bac)
#residuals look fine. 
shapiro.test(mod$residuals)
hist(mod$residuals)
summary(mod)
#Summary of model shows site and timepoint significant predictors
#Interaction of site and timepoint

#test of heterogentity 
leveneTest(Shannon ~ Site, data = richness_BB16S_bac)
leveneTest(Shannon ~ Treatment, data = richness_BB16S_bac)
leveneTest(Shannon ~ Timepoint, data = richness_BB16S_bac)

#significant Levene's test for Site. We will run an additional Kruskal Wallis for site. 
kruskal.test(Shannon ~ Site, data = richness_BB16S_bac)

#What about pairwise differences in site
TukeyHSD(mod, "Site")
plot(TukeyHSD(mod, "Site"), las=1, cex.axis = 0.4)

#Several significant pairwise differences to accompany the significant global test for interaction term. We can
TukeyHSD(mod, "Site:Timepoint")
plot(TukeyHSD(mod, "Site:Timepoint"), las=1, cex.axis = 0.4)

Tests of Chao 1

#ANOVA and Tukey's PW comparisons
#Run a full 3 way ANOVA with interactions
mod<-aov(Chao1 ~ Site * Treatment * Timepoint, data = richness_BB16S_bac)
#residuals look fine. 
shapiro.test(mod$residuals)
hist(mod$residuals)
summary(mod)

leveneTest(Chao1 ~ Site, data = richness_BB16S_bac)
leveneTest(Chao1 ~ Treatment, data = richness_BB16S_bac)
leveneTest(Chao1 ~ Timepoint, data = richness_BB16S_bac)

#Site and the interactions with site are significant. 

TukeyHSD(mod, "Site")
plot(TukeyHSD(mod, "Site"), las=1, cex.axis = 0.4)

Tests of Observed - unique taxa in a sample.

#ANOVA and Tukey's PW comparisons
#Run a full 3 way ANOVA with interactions
mod<-aov(Observed ~ Site * Treatment * Timepoint, data = richness_BB16S_bac)
#residuals look fine. 
shapiro.test(mod$residuals)
summary(mod)
hist(mod$residuals)

leveneTest(Observed ~ Site, data = richness_BB16S_bac)
leveneTest(Observed ~ Treatment, data = richness_BB16S_bac)
leveneTest(Observed ~ Timepoint, data = richness_BB16S_bac)

#Summary of model shows site and timepoint to be a significant predictor as well as the site interaction terms.
TukeyHSD(mod, "Site")
plot(TukeyHSD(mod, "Site"), las=1, cex.axis = 0.4)

#Beta diversity Full model with no strata argument

#extract data from phyloseq object
tab_adonis<-data.frame(otu_table(bb16S_bac_ML))
sd_adonis<-data.frame(sample_data(bb16S_bac_ML))
sd_adonis$Treatment <-as.factor(sd_adonis$Treatment)
sd_adonis$Site <-as.factor(sd_adonis$Site)
dist_w<-phyloseq::distance(bb16S_bac_ML, method =  "wunifrac")
dist_uw<-phyloseq::distance(bb16S_bac_ML, method =  "unifrac")

#full model with no strata arguemnt. Essentially the effect of each main effect without controlling for any potential groups. 

#weighted
adonis2(dist_w ~  Site * Treatment * Timepoint, data= sd_adonis, permutations = 10000)
#unweighted
adonis2(dist_uw ~  Site * Treatment * Timepoint, data= sd_adonis, permutations = 10000)

#with strata argument
adonis2(dist_w ~  Site * Treatment * Timepoint, data= sd_adonis, strata = sd_adonis$Site, permutations = 10000)
adonis2(dist_uw ~  Site * Treatment * Timepoint, data= sd_adonis, strata = sd_adonis$Site, permutations = 10000)


##PW adonis w/ strata = site
pairwise.adonis2(dist_uw ~ Treatment + Timepoint, strata = "Site",  data= sd_adonis, nperm = 10000)

Ordinations

#weighted Unifrac
ord_nmds<-ordinate(physeq = bb16S_bac_ML, method = "NMDS", distance = "wunifrac")
ord_nmds$stress
plot_ordination(physeq = bb16S_bac_ML, ordination = ord_nmds, type = "Samples", color = "Site", shape = "Treatment" ) + theme_classic() + ggtitle("Weighted Unifrac NMDS") + theme(plot.title = element_text(hjust = 0.5, size = 18)) + geom_point(size = 3) + theme(axis.text = element_text(size = 14), axis.title=element_text(size=14))
#ggsave("../../Writing/Final_Figs/Weighted_Unifrac_NMDS.eps")
#strongest grouping by Site. Time point also has a smaller effect, but these groupings are more clear than in the CAP ordination. There are several outlier samples that ordinate far to the right on the NMDS1 axis. 

#unweighted Unifrac
ord_nmds<-ordinate(physeq = bb16S_bac_ML, method = "NMDS", distance = "unifrac")
ord_nmds$stress
plot_ordination(physeq = bb16S_bac_ML, ordination = ord_nmds, type = "Samples", color = "Site", shape = "Treatment" ) + theme_classic() + ggtitle("Unweighted Unifrac NMDS") + theme(plot.title = element_text(hjust = 0.5, size = 18)) + geom_point(size = 3) + theme(axis.text = element_text(size = 14), axis.title=element_text(size=14))
#ggsave("../../Writing/Final_Figs/Unweighted_Unifrac_NMDS.eps")

Sampling time 2 (post treatment)

#post-treatment samples
bb16S_bac_ML_B<-subset_samples(bb16S_bac_ML, Timepoint == "B")

#extract data from phyloseq object for post-treatment
sd_adonis<-data.frame(sample_data(bb16S_bac_ML_B))
sd_adonis$Treatment <-as.factor(sd_adonis$Treatment)
sd_adonis$Site <-as.factor(sd_adonis$Site)
dist_wb<-phyloseq::distance(bb16S_bac_ML_B, method =  "wunifrac")
dist_uwb<-phyloseq::distance(bb16S_bac_ML_B, method =  "unifrac")


#with strata argument
adonis2(dist_wb ~  Site * Treatment, data= sd_adonis, strata = sd_adonis$Site, permutations = 10000)
adonis2(dist_uwb ~  Site * Treatment, data= sd_adonis, strata = sd_adonis$Site, permutations = 10000)

pairwise.adonis2(dist_uwb ~Treatment, strata = "Site", data = sd_adonis, nperm = 10000)

#ordinations
#weighted Unifrac
ord_nmds<-ordinate(physeq = bb16S_bac_ML_B, method = "NMDS", distance = "wunifrac")
ord_nmds$stress
plot_ordination(physeq = bb16S_bac_ML_B, ordination = ord_nmds, type = "Samples", color = "Site", shape = "Treatment" ) + theme_classic() + ggtitle("Weighted Unifrac NMDS") + theme(plot.title = element_text(hjust = 0.5, size = 18)) + geom_point(size = 3)
ggsave("../../Writing/Final_Figs/Time2_Weighted_Unifrac_NMDS.pdf")
#strongest grouping by Site. Time point also has a smaller effect, but these groupings are more clear than in the CAP ordination. There are several outlier samples that ordinate far to the right on the NMDS1 axis. 

#unweighted Unifrac
ord_nmds<-ordinate(physeq = bb16S_bac_ML_B, method = "NMDS", distance = "unifrac")
ord_nmds$stress
plot_ordination(physeq = bb16S_bac_ML_B, ordination = ord_nmds, type = "Samples", color = "Site", shape = "Treatment" ) + theme_classic() + ggtitle("Unweighted Unifrac NMDS") + theme(plot.title = element_text(hjust = 0.5, size = 18)) + geom_point(size = 3)
ggsave("../../Writing/Final_Figs/Time2_Unweighted_Unifrac_NMDS.pdf")

create taxa bar plot to display makeup

sort(table(data.frame(tax_table(bb16S_bac_rarefy))$Phylum), decreasing = TRUE)
length(unique(data.frame(tax_table(bb16S_bac_rarefy))$Phylum))

sort(table(data.frame(tax_table(bb16S_bac_rarefy))$Order), decreasing = TRUE)
length(unique(data.frame(tax_table(bb16S_bac_rarefy))$Order))

sort(table(data.frame(tax_table(bb16S_bac_rarefy))$Family), decreasing = TRUE)
length(unique(data.frame(tax_table(bb16S_bac_rarefy))$Family))

plot_bar(subset_samples(bb16S_bac_rarefy, Timepoint == "A"), fill = "Phylum")  + scale_fill_manual(values = col_vector) + geom_bar(stat = "identity")+ coord_flip()
ggsave("../../Writing/Final_Figs/taxabarA.pdf")

plot_bar(subset_samples(bb16S_bac_rarefy, Timepoint == "B"), fill = "Phylum")  + scale_fill_manual(values = col_vector) + geom_bar(stat = "identity")+ coord_flip()
ggsave("../../Writing/Final_Figs/taxabarB.pdf")

#Core analysis

#create super long color vector
col_vector <- c("#000000", "#FFFF00", "#1CE6FF", "#FF34FF", "#FF4A46", "#008941", "#006FA6", "#A30059",
        "#FFDBE5", "#7A4900", "#0000A6", "#63FFAC", "#B79762", "#004D43", "#8FB0FF", "#997D87",
        "#5A0007", "#809693", "#FEFFE6", "#1B4400", "#4FC601", "#3B5DFF", "#4A3B53", "#FF2F80",
        "#61615A", "#BA0900", "#6B7900", "#00C2A0", "#FFAA92", "#FF90C9", "#B903AA", "#D16100",
        "#DDEFFF", "#000035", "#7B4F4B", "#A1C299", "#300018", "#0AA6D8", "#013349", "#00846F",
        "#372101", "#FFB500", "#C2FFED", "#A079BF", "#CC0744", "#C0B9B2", "#C2FF99", "#001E09",
        "#00489C", "#6F0062", "#0CBD66", "#EEC3FF", "#456D75", "#B77B68", "#7A87A1", "#788D66",
        "#885578", "#FAD09F", "#FF8A9A", "#D157A0", "#BEC459", "#456648", "#0086ED", "#886F4C",
        
        "#34362D", "#B4A8BD", "#00A6AA", "#452C2C", "#636375", "#A3C8C9", "#FF913F", "#938A81",
        "#575329", "#00FECF", "#B05B6F", "#8CD0FF", "#3B9700", "#04F757", "#C8A1A1", "#1E6E00",
        "#7900D7", "#A77500", "#6367A9", "#A05837", "#6B002C", "#772600", "#D790FF", "#9B9700",
        "#549E79", "#FFF69F", "#201625", "#72418F", "#BC23FF", "#99ADC0", "#3A2465", "#922329",
        "#5B4534", "#FDE8DC", "#404E55", "#0089A3", "#CB7E98", "#A4E804", "#324E72", "#6A3A4C",
        "#83AB58", "#001C1E", "#D1F7CE", "#004B28", "#C8D0F6", "#A3A489", "#806C66", "#222800",
        "#BF5650", "#E83000", "#66796D", "#DA007C", "#FF1A59", "#8ADBB4", "#1E0200", "#5B4E51",
        "#C895C5", "#320033", "#FF6832", "#66E1D3", "#CFCDAC", "#D0AC94", "#7ED379", "#012C58",
        
        "#7A7BFF", "#D68E01", "#353339", "#78AFA1", "#FEB2C6", "#75797C", "#837393", "#943A4D",
        "#B5F4FF", "#D2DCD5", "#9556BD", "#6A714A", "#001325", "#02525F", "#0AA3F7", "#E98176",
        "#DBD5DD", "#5EBCD1", "#3D4F44", "#7E6405", "#02684E", "#962B75", "#8D8546", "#9695C5",
        "#E773CE", "#D86A78", "#3E89BE", "#CA834E", "#518A87", "#5B113C", "#55813B", "#E704C4",
        "#00005F", "#A97399", "#4B8160", "#59738A", "#FF5DA7", "#F7C9BF", "#643127", "#513A01",
        "#6B94AA", "#51A058", "#A45B02", "#1D1702", "#E20027", "#E7AB63", "#4C6001", "#9C6966",
        "#64547B", "#97979E", "#006A66", "#391406", "#F4D749", "#0045D2", "#006C31", "#DDB6D0",
        "#7C6571", "#9FB2A4", "#00D891", "#15A08A", "#BC65E9", "#FFFFFE", "#C6DC99", "#203B3C",
        "#671190", "#6B3A64", "#F5E1FF", "#FFA0F2", "#CCAA35", "#374527", "#8BB400", "#797868",
        "#C6005A", "#3B000A", "#C86240", "#29607C", "#402334", "#7D5A44", "#CCB87C", "#B88183",
        "#AA5199", "#B5D6C3", "#A38469", "#9F94F0", "#A74571", "#B894A6", "#71BB8C", "#00B433",
        "#789EC9", "#6D80BA", "#953F00", "#5EFF03", "#E4FFFC", "#1BE177", "#BCB1E5", "#76912F",
        "#003109", "#0060CD", "#D20096", "#895563", "#29201D", "#5B3213", "#A76F42", "#89412E",
        "#1A3A2A", "#494B5A", "#A88C85", "#F4ABAA", "#A3F3AB", "#00C6C8", "#EA8B66", "#958A9F",
        "#BDC9D2", "#9FA064", "#BE4700", "#658188", "#83A485", "#453C23", "#47675D", "#3A3F00",
        "#061203", "#DFFB71", "#868E7E", "#98D058", "#6C8F7D", "#D7BFC2", "#3C3E6E", "#D83D66",
        "#2F5D9B", "#6C5E46", "#D25B88", "#5B656C", "#00B57F", "#545C46", "#866097", "#365D25",
        "#252F99", "#00CCFF", "#674E60", "#FC009C", "#92896B")

calculate core at genus and asv level using rarefied data

Core Genus

#core asv adonis
tab_adonis<-data.frame(otu_table(core_genus))
sd_adonis<-data.frame(sample_data(core_genus))
sd_adonis$Treatment <-as.factor(sd_adonis$Treatment)
sd_adonis$Site <-as.factor(sd_adonis$Site)
dist_core_w_genus<-phyloseq::distance(core_genus, method =  "wunifrac")
dist_core_uw_genus<-phyloseq::distance(core_genus, method =  "unifrac")

#weighted unifrac
adonis2(dist_core_w_genus ~  Site * Treatment * Timepoint, data= sd_adonis, strata = sd_adonis$Site, permutations = 1000)
#unweighted unifrac
adonis2(dist_core_uw_genus ~  Site * Treatment * Timepoint, data= sd_adonis, strata = sd_adonis$Site, permutations = 1000)

#ordinations
#weighted Unifrac
ord_nmds<-ordinate(physeq = core_genus, method = "NMDS", distance = "wunifrac")
ord_nmds$stress
plot_ordination(physeq = core_genus, ordination = ord_nmds, type = "Samples", color = "Site", shape = "Treatment" ) + theme_classic() + ggtitle("Core Genera Weighted Unifrac NMDS") + theme(plot.title = element_text(hjust = 0.5, size = 18)) + geom_point(size = 3)  + theme(axis.text = element_text(size = 14), axis.title=element_text(size=14))
ggsave("../../Writing/Final_Figs/CoreGenus_Weighted_Unifrac_NMDS.eps")
#strongest grouping by Site. Time point also has a smaller effect, but these groupings are more clear than in the CAP ordination. There are several outlier samples that ordinate far to the right on the NMDS1 axis. 

#unweighted Unifrac
ord_nmds<-ordinate(physeq = core_genus, method = "NMDS", distance = "unifrac")
ord_nmds$stress
plot_ordination(physeq = core_genus, ordination = ord_nmds, type = "Samples", color = "Site", shape = "Treatment" ) + theme_classic() + ggtitle("Core Genera Unweighted Unifrac NMDS") + theme(plot.title = element_text(hjust = 0.5, size = 18)) + geom_point(size = 3) + theme(axis.text = element_text(size = 14), axis.title=element_text(size=14))
ggsave("../../Writing/Final_Figs/CoreGenus_Unweighted_Unifrac_NMDS.eps")

Visualization of core microbiome using barplot

#pull out core genera names
core_genus_names<-rownames(tax_table(core_genus))
genus_OTU_table<-data.frame(t(otu_table(bb16S_bac_rarefy_tax_genus)))
#check orientation of otu table
rownames(genus_OTU_table)
#add columns corresponding to core status, % sites, % reads
genus_OTU_table$Core <- NA 
for(i in 1:nrow(genus_OTU_table)){
  if(rownames(genus_OTU_table)[i] %in% core_genus_names == T){
  genus_OTU_table[i,]$Core <- "Core"
  } else{
  genus_OTU_table[i,]$Core <- "NonCore"
  }
}

#Prevalence column
genus_OTU_table$prev <- NA
for(i in 1:nrow(genus_OTU_table)){
  genus_OTU_table[i,]$prev<-table(genus_OTU_table[i,(1:69)] > 0)["TRUE"]/ 69}

#proportional abundance
genus_OTU_table$abund <- NA
  for(i in 1:nrow(genus_OTU_table)){
  genus_OTU_table[i,]$abund <- rowSums(genus_OTU_table[i,c(1:69)])/sum(colSums(genus_OTU_table[,c(1:69)]))
  }

dat<-genus_OTU_table[,c(70:72)]
dat$otu <- as.factor(rownames(dat))
dat$otu <- reorder(dat$otu, -dat$prev )
dat$Phylum <- as.factor(data.frame(tax_table(bb16S_bac_rarefy_tax_genus))$Phylum)
dat$number <- 1:nrow(dat)

ggplot(dat, aes(x = otu, y=prev, fill = Phylum)) +
  geom_bar(stat="identity") + scale_fill_manual(values=col_vector) + theme_classic() +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) + ylab("Prevalance") + theme(legend.position = "none")
#ggsave("../../Writing/Final_Figs/core_prev_owlegend.eps")

ggplot(dat, aes(x = otu, y=abund,  fill = Phylum)) +
  geom_bar(stat="identity") + scale_fill_manual(values=col_vector) + theme_classic()  + scale_y_reverse() + theme(axis.text.x = element_text(size = 4 , angle = 45, hjust = 1))
#ggsave("../../Writing/Final_Figs/core_abund.eps")

ggplot(dat, aes(x = otu, y=abund,  fill = Phylum)) +
  geom_bar(stat="identity") + scale_fill_manual(values=col_vector) + theme_classic()  + scale_y_reverse() + theme(axis.text.x = element_text(size = 4 , angle = 45, hjust = 1)) + theme(legend.position = "none") + ylab("Proportional Abundance") + xlab("OTU ID")
#ggsave("../../Writing/Final_Figs/core_abund_wolegend.eps")


#color by core inclusion

ggplot(dat, aes(x = otu, y=prev, fill = Core)) +
  geom_bar(stat="identity") + scale_fill_manual(values=c("#88CCEE", "#999933")) + 
  theme_classic() +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) + ylab("Prevalance") + theme(legend.position = "none")  + theme(axis.text.y = element_text(size=12), axis.title.y = element_text(size = 16))
ggsave("../../Writing/Final_Figs/core_prev_owlegend_Core.eps")

ggplot(dat, aes(x = otu, y=abund,  fill = Core)) +
  geom_bar(stat="identity") + scale_fill_manual(values=c("#88CCEE", "#999933")) + 
  theme_classic()  + scale_y_reverse() + theme(axis.text.x = element_text(size = 4 , angle = 45, hjust = 1, face = "bold"))+ theme( axis.text.x=element_blank(), axis.ticks.x=element_blank())  + theme(axis.text.y = element_text(size=12), axis.title.y = element_text(size = 16))
ggsave("../../Writing/Final_Figs/core_abund_Core.eps")

ggplot(dat, aes(x = otu, y=abund,  fill = Core)) +
  geom_bar(stat="identity") +  scale_fill_manual(values=c("#88CCEE", "#999933")) +
  theme_classic()  + scale_y_reverse() +# theme(axis.text.x = element_text(size = 4 , angle = 45, hjust = 1)) + 
theme(legend.position = "none") + ylab("Proportional Abundance") + xlab("OTU") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank()) + theme(axis.text.y = element_text(size=12), axis.title.y = element_text(size = 16))
ggsave("../../Writing/Final_Figs/core_abund_wolegend_Core.eps")

Acidophilic bacteria

analysis of lactic acid bacteria and acetic acid bacteria at time pont 2 (following application)

bb16S_bac_ML_B<-subset_samples(bb16S_bac_ML, Timepoint == "B")
#ACETO
#36 taxa
acetobacter <- subset_taxa(bb16S_bac_ML_B, Family =="Acetobacteraceae")
acetobacter_glomed<- tax_glom(acetobacter, taxrank = "Family")

Warning in FUN(X[[i]], ...) :
  merge_taxa attempted to reduce tree to 1 or fewer tips.
 tree replaced with NULL.

#anova for all taxa combined
dat<-data.frame(otu_table(acetobacter_glomed))
samp_d<-data.frame(sample_data(acetobacter_glomed))
merged_dat<-cbind(dat, samp_d)
#check to make sure md and otu table are in same order. All should be true. 
table(merged_dat$X == rownames(merged_dat))

TRUE 
  36 

##mod<-aov(OTU_40 ~ Treatment, data = merged_dat)
#summary(mod)
#plot(residuals(mod))
#boxplot(merged_dat$OTU_40 ~ merged_dat$Treatment)

#baruta for individual taxa
aceto_df <- cbind(data.frame(otu_table(acetobacter)), data.frame(sample_data(acetobacter)))
table(aceto_df$X == rownames(aceto_df))

TRUE 
  36 

boxplot(aceto_df$OTU_82 ~ aceto_df$Treatment)

View(tax_table(acetobacter))
#otu 82 ID'ed by baruta. More abundant in ORG. Taxonomy for this taxa "Acetobacteraceae, Bombella"

for(i in 1:length(names(data.frame(otu_table(acetobacter))))){
  boxplot(unlist(aceto_df[i]) ~ aceto_df$Treatment)
}

#very zero inflated so we will not use anova. 
#for(i in 1:length(names(data.frame(otu_table(acetobacter))))){
#  mod<-aov(unlist(aceto_df[i]) ~ aceto_df$Treatment)
#  summary(mod)
#  print(summary(mod))
#  print(names(aceto_df[i]))
#}

for(i in 1:length(names(data.frame(otu_table(acetobacter))))){
  mod<-kruskal.test(unlist(aceto_df[i]) ~ aceto_df$Treatment)
  print(mod)
  print(names(aceto_df[i]))
}

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = 0.79956, df = 2, p-value = 0.6705

[1] "OTU_40"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = 1.1129, df = 2, p-value = 0.5732

[1] "OTU_77"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = 6.2678, df = 2, p-value = 0.04355

[1] "OTU_82"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = 3.9036, df = 2, p-value = 0.142

[1] "OTU_101"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = 2.3717, df = 2, p-value = 0.3055

[1] "OTU_122"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = 3.7257, df = 2, p-value = 0.1552

[1] "OTU_152"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = 0.056997, df = 2, p-value = 0.9719

[1] "OTU_181"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = 0.75967, df = 2, p-value = 0.684

[1] "OTU_282"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = 1.031, df = 2, p-value = 0.5972

[1] "OTU_694"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = NaN, df = 2, p-value = NA

[1] "OTU_714"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = 2, df = 2, p-value = 0.3679

[1] "OTU_1129"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = 2, df = 2, p-value = 0.3679

[1] "OTU_1246"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = 2, df = 2, p-value = 0.3679

[1] "OTU_1289"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = 2, df = 2, p-value = 0.3679

[1] "OTU_1322"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = NaN, df = 2, p-value = NA

[1] "OTU_1346"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = NaN, df = 2, p-value = NA

[1] "OTU_1374"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = 2, df = 2, p-value = 0.3679

[1] "OTU_1446"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = 2, df = 2, p-value = 0.3679

[1] "OTU_1522"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = NaN, df = 2, p-value = NA

[1] "OTU_1608"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = 2, df = 2, p-value = 0.3679

[1] "OTU_1628"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = NaN, df = 2, p-value = NA

[1] "OTU_1647"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = NaN, df = 2, p-value = NA

[1] "OTU_1657"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = 2, df = 2, p-value = 0.3679

[1] "OTU_1702"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = 2, df = 2, p-value = 0.3679

[1] "OTU_1709"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = 2, df = 2, p-value = 0.3679

[1] "OTU_1755"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = 2, df = 2, p-value = 0.3679

[1] "OTU_1756"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = 2, df = 2, p-value = 0.3679

[1] "OTU_1809"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = NaN, df = 2, p-value = NA

[1] "OTU_1816"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = 2, df = 2, p-value = 0.3679

[1] "OTU_1919"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = NaN, df = 2, p-value = NA

[1] "OTU_1928"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = NaN, df = 2, p-value = NA

[1] "OTU_2017"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = 2, df = 2, p-value = 0.3679

[1] "OTU_2037"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = 2, df = 2, p-value = 0.3679

[1] "OTU_2215"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = 2, df = 2, p-value = 0.3679

[1] "OTU_2255"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = NaN, df = 2, p-value = NA

[1] "OTU_2286"

    Kruskal-Wallis rank sum test

data:  unlist(aceto_df[i]) by aceto_df$Treatment
Kruskal-Wallis chi-squared = 2, df = 2, p-value = 0.3679

[1] "OTU_2306"

#OTU 82 ID'ed as significantly different among groups
pairwise.wilcox.test(aceto_df$OTU_82, aceto_df$Treatment)

Warning in wilcox.test.default(xi, xj, paired = paired, ...) :
  cannot compute exact p-value with ties
Warning in wilcox.test.default(xi, xj, paired = paired, ...) :
  cannot compute exact p-value with ties
Warning in wilcox.test.default(xi, xj, paired = paired, ...) :
  cannot compute exact p-value with ties

    Pairwise comparisons using Wilcoxon rank sum test with continuity correction 

data:  aceto_df$OTU_82 and aceto_df$Treatment 

    CF    CON  
CON 0.190 -    
ORG 0.058 0.241

P value adjustment method: holm 

boxplot(aceto_df$OTU_82~ aceto_df$Treatment)

ggplot(aceto_df, aes(x = Treatment, y = OTU_82, color = Treatment)) + geom_boxplot()  + theme_classic()  + theme(axis.text = element_text(size = 14), axis.title=element_text(size=14))
ggsave("/Users/gordoncuster/Desktop/Git_Projects/BeeBread/Writing/Final_Figs/OTU_82_ML.eps")

Saving 7.29 x 4.51 in image

#no significant differences across treatment types

##LACTO

#75 taxa
lactobac <- subset_taxa(bb16S_bac_ML_B, Order =="Lactobacillales")
lactobac_glomed <- tax_glom(lactobac, taxrank = "Order")

dat<-data.frame(otu_table(lactobac_glomed))
samp_d<-data.frame(sample_data(lactobac_glomed))
merged_dat<-cbind(dat, samp_d)
#check to make sure md and otu table are in same order. All should be true. 
table(merged_dat$X == rownames(merged_dat))

lacto_df <- cbind(data.frame(otu_table(lactobac)), data.frame(sample_data(lactobac)))
table(lacto_df$X == rownames(lacto_df))
boxplot(lacto_df$OTU_530 ~ lacto_df$Treatment)
View(tax_table(lactobac))



for(i in 1:length(names(data.frame(otu_table(lactobac))))){
  boxplot(unlist(lacto_df[i]) ~ lacto_df$Treatment)
}

for(i in 1:length(names(data.frame(otu_table(lactobac))))){
   mod<-kruskal.test(unlist(lacto_df[i]) ~ lacto_df$Treatment)
  print(mod)
  print(names(lacto_df[i]))
}

https://stackoverflow.com/questions/61796484/stacked-bar-chart-with-count-above-and-below-x-axis